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PACS 64.70.D- - Solid-liquid transitions 

PACS 64.60. qe - General theory and computer simulations of nucleation 
PACS 64.60. an - Finite-size systems 

Abstract -We report an extensive Monte-Carlo study of the melting of the classical two dimen¬ 
sional Wigner crystal for a system of point particles interacting via the 1/r-Coulomb potential. 
A hexatic phase is found in systems large enough. With the multiple histograms method and the 
finite size scaling theory, we show that the fluid/hexatic phase transition is weakly first order. 
No set of critical exponents, consistent with a Kosterlitz-Thouless transition and the finite size 
scaling analysis for this transition, have been found. 


Like charged particles immersed in a homogeneous neu- 
'] H tralizing background form crystals at low temperature 
•301 ; these crystals are called Wigner crystals after the 
Ch seminal work of E.P. Wigner in 1934 on electrons in metals 
0 11]. Since the original work of Wigner, Coulomb crystals 
i_have been observed in a large variety of systems: in plasma 
physics [4j[5], in colloids science |6|[7], in semiconductors 
[8-10 and in biology [7 . Two dimensional Wigner crys¬ 


tals are observed in complex plasmas 


in electrons 
in laser- 
in inver- 


[— trapped on the surface of liquid Helium [11} - 
l/-) cooled 9 Be + ions confined in Penning traps 15 
cn sion layer of semiconductors at low temperature [16 

Colloids and dusty plasmas are classical systems ; the clas- 
sical regime for electrons and ions is defined, in absence of 
CD magnetic field, when the Fermi energy is small compared 
to interaction energy and temperature ; for surface elec- 
• .trons, it corresponds to low surface density [3,17,18]. 

The ground-state of the classical two dimensional Wigner 
crystal is known to be a triangular lattice [2]|3j[T9}|21| and 
$__i it is worthwhile to outline that the long ranged nature of 
73 the interaction in Coulomb systems does not fulfill the hy¬ 
pothesis of the Mermin theorem on the abscence of long 
ranged cristalline order in two dimensions 


22 


For all experimental systems mentioned above, the study 
of the phase transition between the ordered Wigner crys¬ 
tal and the disordered fluid-like phases has peculiar impor¬ 
tance. Not only for a better understanding of the struc¬ 
tural properties of these systems, but also for the theoret¬ 


ical study of the two dimensional meltings 23 - 40 


In this letter, we report an extensive Monte Carlo study 


of the melting of the classical two dimensional Wigner 
crystal. The system is made of N charged point particles 
confined in a two dimensional plane ; periodic boundary 
conditions in two dimensions are used. The interaction 
energy between a pair of particles is the Coulomb energy 
V[r) = Q 2 /r with Q the charge of particles and r the 
distance in the plane between both particles. The charges 
of particles are neutralized by an uniform background of 
charge density <jq ; electroneutrality of the system reads as 
NQ + <r 0 S = 0 with S the surface of the simulation cell. 
The number density is noted p = N/S and ao = —Qp- 
This system is a One Component Plasma (OCP) confined 
to a plane. 

With these notations, for a given configuration {r,}i<j<jv 
of the charged particles, the total energy of the system is 
given by 
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with S 0 the simulation cell and S n , its periodic image. 
The lattice sums in Eqjl] are computed with the Ewald 


method 41 


The Monte Carlo simulations are done at finite tempera- 
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Fig. 1: Thermodynamical properties of the one component 
plasma confined in a plane (2D), (a) Excess internal energy 
per particle as function of the coupling constant T. The sym¬ 
bols in blue with error bars are averages energies computed in 
Monte-Carlo (MC) simulations. The red line is the internal 
energy of the 2D Wigner crystal in the low temperature limit 
as predicted by the harmonic approximation Eq. ©• (b) De¬ 
viation of the excess free energy per particle from the Wigner 
crystal free energy as function of T ; /3Af(T) is computed with 
Eq-©. Inset : enlargement near the phase transition. 


ture T in the canonical ensemble with a variable shape of 
Sq, but at constant surface area. The trial move for the 
shape of the box with the Ewald method is described in 
ref. 42 ; it is particularly well adapted to study solid-solid 
and solid-liquid transitions and it had been used for the 
study of the crystal phases of Coulomb 42 and Yukawa 


bilayers 43 


The only relevant thermodynamical variable in the classi¬ 
cal regime is the coupling constant T = Q 2 ^Wp/ksT. 

The numerical simulations are performed on systems with 
N = 1024, 2025, 4096 and 8100 particles ; for each system 
size, we have performed MC simulations for about 20 dif¬ 
ferent values of T ranging from 126 to 214. The analysis of 
data is done with the multiple histogram method (MHM) 
[44]|47j and the finite size scaling theory 47 54]. 


We define a MC-cycle as a trial move of N particles and a 
trial change of the shape of the simulation box. Equilibra¬ 
tions for each coupling constants are done with 7.5 x 10 5 
MC-cycles for systems with N = 1024 and 2 x 10 5 MC- 


cycles for N = 8100. In the Monte Carlo sampling of the 
phase space, after equilibration, the averages are taken 
over 2.5 x 10 6 MC-cycles for systems with N = 1024 
and with 3.0 — 7.5 x 10 5 MC-cycles for N = 8100. For 
systems size N = 2025 and 4096, the numbers of MC- 
cycles used for equilibrations and averages are in the same 
range. With the code and the parameters of the Ewald 
method used in this study, the typical cpu-times to per¬ 
form 2.0 x 10 4 MC-cycles for the system with N = 8100 
is about 48 hours ; this computing time includes the com¬ 
putational effort for the Voronoi constructions 55 done 
every MC-cycle. 

After equilibration, the trajectories for all systems size 
and all coupling constants are saved to permit the imple¬ 
mentation of the multiple histogram reweighting method 
(MHM) [44-47 and the finite size analysis of the phase 


transition 147f |54| ; we define the linear size of systems as 


L = \fN~fp. 


In refs. 56 57 , it is shown that to observe a stable hex- 


atic phase one needs simulations long enough such as the 
phase space of the system is correctly sampled ; this is re¬ 
lated to the critical properties of hexatic phases 26 . For 


Lennard-Jones systems, the stability of the hexatic phase 
is confirmed by longer computations 


58 


An estimate of the sampling of the phase space is given 
by the average length of the trajectory of particles. In 
the Molecular Dynamics reported in refs. 56 57 , we may 


estimate the average length of trajectories from the aver¬ 
age velocity and the total duration of the computations. 
From the data reported in these works, we found that the 
lengths of trajectories are 25-70 L. In the MC reported 
in the present work, with an amplitude of the trial moves 
as 0.004 — 0.006 L , the average length of trajectories are 
about 550 L for equilibration and 2000 L for averages in 
systems with N = 1024, and about 130 L for equilibration 
and 440 L for averages in systems with N = 8100. 

In the following, we note < . > the averages obtained 
with the MHM and by < . >mc the averages computed 
in Monte Carlo simulations. 

Figure [l] gives the thermodynamical properties of the OCP 
confined in a plane. The excess internal energy per parti¬ 
cles is computed as 


<E>=--^]nZ N (T) 


( 2 ) 


where lnZjv(T) is the partition function computed with 
MHM. The pressure is related to < E > and the contri¬ 
bution to energy due to the neutralizing background (see 
e.g. refs. [59,60].) 

For large value of T, the low temperature limit, the sys¬ 
tem form a triangular Wigner crystal, the harmonic ap¬ 
proximation allows to represent the energy in this limit as 


61 62 


<E> 1 

NQ 2 ^/p - eo + Co r 


(3) 


Letting eo as a free parameter, we find that it does not 
depend on the number of particles and we have eo = 


P-2 
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— 1.961(1), this value is very close to the known Madelung 
constant (cm = —1.960515788...) of the triangular lattice 

H- 

If we fix eo = Cm , then we obtain, for all system sizes, 
Co = 1.845(5). In Figjl] (a), we represent Eq.([3]) with 
eo = Cm and Co = 1.845(5) as a thick solid red line. 
Eq.([2j), with the MHM, permits to obtain the ground state 
energy of a model system to a very good accuracy 1 46 ; 
therefore, it provides a reference point for the computa¬ 
tion of the free energy. 

The partition function Z S (T) in the low temperature limit 
is obtained by integration of Eq.([3]) and, with the defini¬ 
tion given in Eq.([2j), we find 

1 In Z s (T) = Z 0 - ' (c M r + Co In T) (4) 
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where Z$ is a constant of integration and the prefactor 
(2p/7rv / 3 ) 1 ^ 2 stems from the definition of T and the geo¬ 


metric features of the triangular lattice 21 61 


The deviation of the excess free energy per particle from 
the Wigner crystal free energy as function of T is repre¬ 
sented in Figjl] (b) ; it is computed as 

/3A/(r) = 1 InZ n (T) - i In Z S (T) (5) 


where ln^jv(r) is the partition function obtained with 
MHM. 

The slope discontinuity of the tangents near the phase 
transition is an indication of a first order phase transition. 
Since the only relevant thermodynamical parameter is F, 
the first order transition is driven by the change of the 


temperature 

48 
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given by 
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(< E 2 > - < E > 2 ) 


( 6 ) 


<E 4 > 


(7) 


3(< E 2 >) 2 

The variations of Vl with F for different system sizes is 


consistent with a first oder phase transition 47 48 51 


more precisely, the very small deviation of Vl from 2/3 
indicate a weak first order transition 47 . The extremums 
of C/ks and Vl give r c (L), the rounding coupling con¬ 
stant due to the finite size of the simulated systems. We 
define the reduced temperature as 


1 1 

f " r c (L) 


( 8 ) 


and in Figj2]we represent the collpase of the observables 
C/ks (a) and Vl (b) as functions of 7 L 2 . 

The scaling with 7 L 2 is excellent, however the amplitude 
of the thermodynamical observables scales with L and not 



Fig. 2: Finite-size scaling of the thermodynamical observables 
for the one component plasma confined in a plane (2D). All 
data represented are computed with the multiple histogram 
method (MHM). (a) Scaling of the specific-heat C/ks ; inset 
: specific-heat as function of coupling constant T for different 
sizes, (b) Scaling of the fourth-order cumulants of the energy 
; inset : cumulants as function of T. 


with L 2 as excepted for a strong first order phase transi¬ 
tion 48|49,51 . We interpret this anomalous scaling of the 


thermodynamical observables as induced by the weakness 
of the first order phase transition between the liquid and 
hexatic phases 50 . A similar deviation to the finite size 


scaling of first order transition is observed in five-state 
Potts model in two dimensions 50 ; the pseudocritical 


behavior stems from a finite, but very large (> L), corre¬ 
lation length (see e.g. Fig. 16 in ref. 50 ). 

The data for N = 8100 deviates significantly from the 
scaling functions. Periodic boundary conditions favor the 
crystal phases. In small systems (N = 1024 and 2025), the 
range of stability of the hexatic phase with short ranged 
positional order and long ranged orientational order is ar¬ 
tificially reduced because of the positional order induced 
by the periodic boundary conditions. 

In multiple histograms method, one is able to compute 
the partition function In Zn(T) for any values of T in the 
range covered by the MC simulations 44 45 47 and to 
interpolate the variations of observables with T. 

In figure g a), we represent the histograms of energy com- 
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Fig. 3: Histograms (a) and trajectories (b) of the energy close 
to the phase transitions for systems with N = 8100. The 
thin black lines are histograms interpolated from the MC- 
histograms with the multiple histogram method. The thick red 
and blue lines are histograms of energies very close to the tran¬ 
sitions, used to estimate F^ and IV The colored histograms 
represented in (a) with thick lines are MC-histograms com¬ 
puted from the energy trajectories shown in (b) with the same 
color. 


puted from the Monte-Carlo trajectories and histograms 
computed with the MHM for N = 8100 ; in figure [3^b), 
we plot the MC-trajectories of the total energy used to 
construct the MC-Histograms in (a). The orange trajec¬ 
tory and histogram are for the liexatic phase very close 
to the hexatic-liquid transition and, the greens are for the 
hexatic phase very close to the solid-hexatic transition. 
For all systems sizes, we may represent the probability dis¬ 
tribution of the total energy as a superposition of gaussian 


distributions as 48 49 51 



Fig. 4: Location of the maximums of the double (2 peaks) 
and triple (3 peaks) peaked histograms compared to excess in¬ 
ternal energies computed with the multiple histogram method 
(IV = 4096,8100). The structure factors S(k) are computed in 
Monte-Carlo simulations in systems with N = 8100 ; they are 
represented as density plots for three values of T. 


perature close to the phase transitions, the gaussians are 
centered at E a (T) = E a + C a 7 - 

The values for E a (T) with a € {s,h,l} (triple peaks) and 
a £ {s, /} (double peaks) are represented on Figj4]with the 
excess internal energy computed with MHM by Eq. (|2|) for 
N = 4096 and 8100. These data permit to locate the tran¬ 
sition coupling constants and the range of stability of the 
hexatic phase, reported in Table [l] 

In Figj4j we report also the structure factor S(k) for the 
systems with N = 8100, for the fluid phase (T = 137.58), 
the hexatic phase (F = 138.96) and the Wigner crystal 
(F = 140.36). In experiments, the structure factors are ob¬ 
tained from diffraction patterns [63]; in simulations, S(k) 
are computed as 

S( k ) = ( p{k)p(-k)) MC ( 10 ) 

with 

N 

p(k) = / ds exp(— ik.r)p(r) = ^ exp(— ik.r n ) ( 11 ) 

71= l 


H n (E; r ) = ^ H a (T) exp (E ) (9) 

with a £ {s, h, 1} (triple peaks) or a £ {s, 1} (double peaks) 
where s stands for the Wigner crystal phase, h the hexatic 
phase and l the disordered liquid phase. 

In the plot of histograms, the various phases are repre¬ 
sented by the peaks and, for a given size, we locate the 
temperature of transition between two phases, 1 /F C (L), 
by requiring that the height of the peaks to be the same. 
Assuming that the specific heats do not vary with the tem- 


The structure factors, shown on Figj4j support the identi¬ 
fication of the three peaks with the three phases : liquid, 
hexatic and solid. 

The bond orientational order parameter </>6 is a suitable 
quantity for the study of the melting of triangular lattices. 
We compute <pQ, the susceptibility \6 and the fourth-order 
cumulant Uq with Voronoi constructions, as described in 
ref. [43] . 

On Fig(5] we represent the data collapses of < >, \6 

and U 6 as functions of the reduced temperature 7 L 2 . The 
finite size scaling depends on the surface of the system (or 
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Table 1: Estimate of the transition coupling constants and la¬ 
tent heat for the melting of the two dimensional Wigner crystal. 
T h (L) is the coupling constant for the fluid/hexatic transition 
; Li/h the latent heat of the fluid/hexatic transition ; AT and 
E S (L) are respectively crude estimates of the stability range of 
the hexatic phase and the hexatic/solid transition. 


N 

I \(L) 

Li/h 

AT 

r s(L) 

1024 

135.0i0.14 

1.2 x 10" 3 

0.16 

135 

2025 

136.7i0.12 

1.2 x 10" 3 

0.19 

137 

4096 

137.5i0.1 

1.4 x 10 " 3 

0.32 

138 

8100 

138.81i0.05 

1.0 x 10" 3 

1.05 

140 


number of particles) that is consistent with a first order 
phase transition. As for the thermodynamical observables, 
data for systems N = 8100 deviate significantly from the 
scaling functions ; we interpret this deviation as result¬ 
ing of an increase of the thermal stability of the hexatic 
phase in larger systems. On the contrary to the scaling of 
the amplitude of C/ks with the system size L observed in 
Figj2j the amplitude of \6 scales with L 2 as predicted for 
a first order transition. 

On FigJHJd), we represent the bond orientational correla¬ 
tion functions g§{r) for the MC-trajectories of Fig{3])b). 

In the KTHNY theory of the two dimensional melting 
23-27 , the finite size scaling of \6 at the liquid-liexatic 
transition (7 > 0 ) is governed by X 6 ~ Q 2 m with £ 6 
the correlation length of the order parameter ; it exhibits 
an essential singularity scaling as exp(b/ / y u ) at the transi¬ 
tion 23 . In the present study, no set of critical exponents, 


consistent with the KTHNY theory or with the XY-model 


23 , have been found to achieve data collapses better, or 


at least as good as, thoses done for a first order phase 
transition 


52-54 


All the results reported in Figs[l][5] support a weak first 
order phase transition for the liquid-hexatic transition in 
the melting of the Wigner crystal. 

Based on the analysis of the triple peaked histograms, we 
report in Table [l] the transition coupling constants for the 
hexatic/solid transition T S (L) and for the weak first order 
liquid/hexatic transition Th{L) 1 found for all systems size. 
The errors bars on values of T h(L) are obtained by neigh¬ 
boring histograms with peaks about the same height. The 
latent heats of the liquid/hexatic transition are computed 
as Li/ h = Ei~E h (see also Fig|4j) ; the uncertainties on the 
numerical values of Lu h are estimated about 10%. The 
small values of the latent heat found for the fluid/hexatic 
transition is another signature of the weakness of the first 
order phase transition ; the same order of magnitude is 


found for hard disks systems 33 


The fluid/hexatic phase transition for the Coulomb system 
studied in the present work is compatible with the grain 
boundary induced melting found in experiments on com¬ 
plex plasmas [5]. This transition has close similarities with 
the melting of hard disk systems 28 -30 and hard spheres 



Fig. 5: Finite-size scaling for the bond orientational order pa¬ 
rameters and derived quantities, (a) Scaling of the bond orien¬ 
tational order parameter cj> q. (b) Scaling of the susceptibility 
Xe- (c) Scaling of the fourth-order cumulant of the order pa¬ 
rameter t/6.(d) Bond orientational correlation functions ge (r) 
for the MC-trajectories of Fig[3|b). 


in slab geometry 32 , but it is clearly different from the 


KTHNY mechanism found for superparamagnetic colloids 


at air-water interface in an external magnetic held 34-38 


As explained before, the stability of the Wigner crystal is 
increased because of the periodic boundary conditions (see 
also ref. |33]) ; nevertheless, the detailled analysis of the 
histograms of energy (Figs {3] and [4]) have permitted to lo¬ 
cate approximatively the hexatic/solid transition (T S (L) 
in Table [T]). Figures [3] and [3] tend to support a first or¬ 
der phase transition for the hexatic/solid transition, in 
agreement with a previous study done by B.K. Clark and 
co-workers on 2D quantum Coulomb systems 18 . How¬ 


ever, with the system sizes considered in the present work, 
we haven’t yet been able to reach a definitive conclusion 
on the nature and order of the hexatic/solid transition. 


This work is supported by Theme 2 - LabEx PALM 
(projet - SLAB - ANR-10-LABX-0039-PALM) and project 
ECOS-Nord C14P01. The author acknowledges the com¬ 
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tion Informatique of Universite Paris-Sud. 
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